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HIGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES 

PROBLEM 

AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG 


Abstract. We develop a high order cut finite element method for the Stokes problem 
based on general inf-sup stable finite element spaces. We focus in particular on composite 
meshes consisting of one mesh that overlaps another. The method is based on a Nitsche 
formulation of the interface condition together with a stabilization term. Starting from 
inf-sup stable spaces on the two meshes, we prove that the resulting composite method is 
indeed inf-sup stable and as a consequence optimal a priori error estimates hold. 


1. Background 

1.1. Introduction. Meshing of complex geometries remains a challenging and time con¬ 
suming task in engineering applications of the finite element method. There is therefore a 
demand for finite element methods based on more flexible mesh constructions. One such 
flexible mesh paradigm is the formulation of finite element methods on composite meshes 
created by letting several meshes overlap each other. This approach enables using combi¬ 
nations of meshes for certain parts of a domain and reuse of meshes for complicated parts 
that may have been difficult and time consuming to construct. 

We consider the case of a composite mesh consisting of one mesh that overlaps another 
mesh which together provide a mesh of the computational domain of interest. This results 
in some elements on one mesh having an intersection with one or several elements on the 
boundary of the other mesh. We denote such elements by cut elements. The interface 
conditions on these cut elements are enforced weakly and consistently using Nitsche’s 
method [18]. 

In this setting ra first developed and analyzed a composite mesh method for elliptic 
second order problem based on Nitsche’s method. In mi, this approach was extended to 
the Stokes problem using suitable stabilization to ensure inf-sup stability of the method. 
Implementation aspects were discussed in detail in [16j . In [Hj a related cut finite element 
method for a Stokes interface problem based on the Pl-iso-P2 element was developed and 
analyzed. 
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Composite mesh techniques using domain decomposition are often called chimera, see 
for example ra, 0 for uses in a finite difference setting or [12] in a finite element set¬ 
ting. The extended finite element method (XFEM) also provides composite mesh handling 
techniques, see for example P ED]. However, the Nitsche method approach using cut 
elements used in this work makes it possible to obtain a consistent and stable formula¬ 
tion while maintaining the conditioning of the algebraic system for both conforming and 
non-conforming high order finite elements. 

In this paper, we consider Stokes flow and device a method based on a stabilized Nitsche 
formulation for enforcement of the interface conditions at the border between the two 
meshes. A specific feature is that we only assume that we have inf-sup stable spaces on 
the two meshes and that the spaces consist of polynomials. We can then show that our 
stabilized Nitsche formulation satisfies an inf-sup condition and as a consequence optimal 
order also a priori error estimates hold. We emphasize that the spaces are arbitrary and 
can be different on the two meshes, in particular, continuous or discontinuous pressure 
spaces as well as higher order spaces can be used. We present extensive numerical results 
for higher order Taylor-Hood elements in two and three spatial dimensions that confirm 
our theoretical results. 

The outline of the paper is as follows: First we review the Stokes problem. Then the 
finite element method is presented by first defining the composite mesh and introducing 
finite element spaces. The method is then analyzed where the inf-sup condition is the main 
result. Finally we present the numerical results and the conclusions. 

1.2. The Stokes Problem. In this section, we review the Stokes problem and state its 
standard weak formulation. We also introduce some basic notation. 

1.2.1. Strong form. Let 12 be a polygonal domain in with boundary <912. The Stokes 
problem takes the form: Find the velocity u : 12 —)■ M. d and pressure p : 12 —)■ R such that 


(1.1) 

—A u + Vp = f 

in 12, 

(1.2) 

div u = 0 

in 12, 

(1,3) 

u = 0 

on c212, 


where f : 12 — > M. d is a given right-hand side. 

1.2.2. Weak form. As usual, let LP(12) denote the standard Sobolev space of order s > 0 
on 12 with norm denoted by || • ||m*(n) and semi-norm denoted by | • | h« (O)• Let L 2 (12) 
denote the L 2 -norm on 12 with norm denoted by || ■ ||q. The corresponding inner products 
are labeled accordingly. 

Introducing the spaces 

(1.4) V= 

(1.5) Q = {q e L 2 (12) : f qdx = 0}, 

J ci 
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with norms ||Di;||q = ||u g) V||n and \\q\\n, the weak form of (1.1) and (1.2) reads: Find 
(u,p) G V x Q such that 

(1.6) a(u,v) + b(u,q) + b(p, v) — l(v) V(v,q)eVxQ, 
where the forms are defined by 

(1.7) a(u,v) = (Du, Dv)q, 

(1.8) b(u,q) = -(di vu,q) n = (u,Vq)n, 

(1-9) l{v) = {f,v) a . 


Remark 1. We obtain the variational problem (1.6) by formally multiplying (1.1) by a 
test function v and (1.2) by a test function —q. 


It is then possible to show that the inf-sup condition 

( 1 . 10 ) 


n n ^ b(v,q) (div v, q) 

l|g||n<sup =sup VqeQ 

vev\\Dv\\n V& v \\Dv n 


holds, from which it follows that there exists a unique solution to (1.6). See [6] for further 
details. 


2. Methods 

2.1. The Composite Mesh. We here present the concepts and notation of the domains 

and meshes used. The main idea is to introduce a background domain which is partially 
overlapped by another domain (the overlapping domain). For each of these domains, we 
mimic the setup of a traditional finite element method in the sense that each domain is 
equipped with a traditional finite element mesh. The two meshes are completely unrelated. 
In particular, the interface between the two meshes is determined by the overlapping do¬ 
main and is not required to match or align with the triangulation of the background 

domain. 

2.1.1. The composite domain. Let the predomains Oj C fl, i = 0,1, be polygonal subdo¬ 
mains of fI in M. d such that Oo U Hi = fl; see Figure [l} Consider the partition 

(2.1) fl = fl 0 Ufli, 

( 2 . 2 ) n Q = n\n u 

(2.3) fli = fii, 

and let T = dfl\\dfl be the interface between the overlapping domain O] and the underlying 
domain fl 0 ; see Figure [lj We make the basic assumption that each flj, i = 0,1, has a 
nonempty interior. We note that implies that there exists a nonempty open set U G fl 
such that T fl U 0. (The set U plays an important role in the proof of Lemma [7] below.) 
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FIGURE 1. The domains fl* and the subdomains (all shaded) sharing the 
interface F. 



Figure 2. The domains f(shaded). 


2.1.2. The composite mesh. For z = 0,1, let be a quasi-uniform mesh on Q, with mesh 
parameter h E (0, h] and let 

(2.4) F h , = {K E t Ki : K n ± 0} 

be the submesh consisting of elements that intersect f Ip, see Figure [3} Note that Kh ,o 
includes elements that partially intersect Q]. We also introduce the notation 

(2.5) n h) i= U K. 

K&K Ki 

Note that Oi = f4,i and h2 0 C see Figure [2j 

We obtain a partition of fl by intersecting the elements with the subdomains: 

i i 

( 2 . 6 ) U /C M n Qi = U {K n Qi-.Ke ic h}i }. 

i=0 i=0 

See also Figure [3] 


2.2. Finite element formulation. In this section, we present the finite element method 
for approximating the weak form (1.6). Some notation will be introduced, but the main 
idea is to assume we have inf-sup stable spaces in each of the subdomains away from the 
interface. Then we are able to formulate a method similar to [TO] and d 


2.2.1. Finite element spaces. For each of the predomains fb with corresponding family of 
meshes A Zh,i we consider velocity and pressure finite element spaces V^i x Qh,i- The spaces 
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FIGURE 3. The meshes /C h,i and /C h,i of the corresponding domains fh and 
f 1h,i- Note that T is not aligned with /C^q. 


do not contain boundary conditions since these will be enforced by the finite element 
formulation. We define 

(2.7) Vh,i X Qh,i ^ 

where i = 0,1 and define 

i 

(2.8) Vh x Q h = V h , x Qh,i- 

i =0 

Note that since the domains Qh,i overlap each other, Vh x Qh is to be understood as a 
collection of function spaces on the overlapping patches i — 0,1. We now make the 
following fundamental assumptions on these spaces: 

Assumption A (Piecewise polynomial spaces). The finite element spaces Vh and Qh 
consist of piecewise polynomials of uniformly bounded degree k and l, respectively. 


Assumption B (Inf-sup stability). The finite element spaces are inf-sup stable restricted 
to a domain bounded away from the interface. More precisely, we assume that for i — 0,1 
and h G (0, h\ there is a domain oj h ^ C such that: 

(a) The set uih,i is a union of elements in JCh,ii see Figure^ 

(b) The inf-sup condition 


(2.9) 


V M (p)IU,, < sup 

vew h 


(divu,p) 


u h ,i 


\Dv\ 


Wh,i 


holds, where X Uh i (p) is the average of p over u>h,i and Wh,% is the subspace of Vh,i defined 
by 


(2.10) W h)0 = {v G Vh,o : v = 0 on Q h ,o \ u h , o}, 

(2.11) W M = V h>1 . 


(c) The set Uh ,o is close to flo in the sense that 

(2.12) Qh,o \wh,o C [/^(r), 8 ~ h , 

where Us( T) = {x G : |p(cc)| < 5} is the tubular neighborhood ofT with thickness 8. 
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^h, 0 ^h, 1 ^h, 1 

Figure 4. The domains c Oh,i (shaded) where inf-sup stability is assumed. 

Note that T is outside 0 - 

Remark 2. The assumptions presented ensure that the polynomial spaces are such that 
certain inverse inequalities hold. More generally, inverse inequalities hold if there is a finite 
set of finite dimensional reference spaces used to construct the element spaces. The use of 
the interpolant in the proof of Lemma [?| could alternatively be handled using an abstract 
approximation property assumption. 


2.2.2. Finite element method. We consider the finite element method: Find (uh,Ph) £ 
V h x Q h such that 

(2.13) A h ((u h ,p h ), (v,q)) = l h (v) V (v, q) G V h x Q h , 


where the forms are defined by 


(2.14) 

A h((u,p),(v,q)) 

= a h (u , v) + b h (u , q ) + b h (v,p) + d h ((u,p), (v , q)), 

(2.15) 

a h (u,v) 

= a h:N (u, v ) + a h , 0 (u , v), 

(2.16) 

a h>N (u,v) 

= (Du,Dv) no + ( Du,Dv) Ql 



- ({(Du) ■ n), [v]) r - ([w], ((Dv) ■ n)) r 



+ (3h~ l ([u], [w]) r , 

(2.17) 

a ht0 {u,v) 

= ( [Du], [Dv]) n M nn!, 

(2.18) 

h(u,q) 

= —(div w, q)n 0 - (di \u,q) Ql + ([n ■ u ], (q)) r , 

(2.19) 

dh((u,p), (v,q)) 

= h 2 (Au - Vp, Av + V?)n M \ Wh>0 , 

(2.20) 

lh(v) 

= (/, v) n - h 2 (f , Av + Vg)n h>0 \ WM - 


Here, n is the unit normal to F exterior to Hi, [u] = V\ — v 0 is the jump at the interface 
T and (■ v) = (vo + Vi)/2 is the average at T (although any convex combination is valid 
TO- The parameter ft > 0 is the Nitsche parameter and must be sufficiently large (see 
for example m) and scales as k 2 , where k is the polynomial degree. Furthermore, h is 
the representative mesh size of the quasi-uniform mesh. In a practical implementation, h 
is evaluated as the local element size. 

A comment on the respective terms may be clarifying: ah,N and bh are the standard 
Nitsche formulation of (1.6) and ah,o is a stabilization of the jump of the gradients across 
T (see S3)' The least-squares type term dh stabilizes the method since we do not assume 
inf-sup stability in all of H 0 . 
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By simple inspection, we note that the method is consistent. We conclude by noting 
that the method satisfies the Galerkin orthogonality. 


Proposition 3 (Galerkin orthogonality). Let (u,p) G V x Q be a weak solution to the 
formulation (1.6) and let {uh,Ph) € VhXQh be the solution to the finite element formulation 
(2.13). Then it holds 


( 2 . 21 ) 


A h ((u,p) - (u h ,p h ),(v h ,q h )) = 0 V(v h ,q h ) G V h x Q h . 


Proof. The result follows from [TO] and noting that ah t o{ u > v h) = dh((upp), (v h , q h )) = 0 
for all (v h , q h ) G V h X Q h . □ 


2.2.3. Approximation properties. We assume that there is an interpolation operator 7: 
VfiQh,i) Vh,i, for i — 0,1, where V)(fC [L 2 {yi h ^)] d is a space of sufficient regularity 
to define the interpolant. For Taylor-Hood elements, we take 7to be the Scott-Zhang 
interpolation operator [TO] , and VfiLlhj) — [L 2 (filh,i)] d ■ For other elements we refer to their 
corresponding papers, for example the Crouzeix-Raviart element [8j, the Mini element jT], 
or the overviews in [3] or [4]. 

The full interpolation operator 7Zh : V —> Vh can now be defined by the use of a linear 
extension operator E : [H s (uh,o)] d —» ]/T s (r2/i )0 )] d , s > 0, such that (Ev)| Wh 0 = v and 

( 2 - 22 ) ll Ev llH«(n h , 0 ) < IMIh-^o)- 

Now, 7T/ t :V^V h is defined by 

(2.23) 7T h V = TT h fiEVo © 7T h) iVi. 

A similar argument can be made to define the pressure interpolation operator 7 Th : Q ^ Qh- 
Furthermore, we assume that the following standard interpolation estimate holds: 

(2.24) \\v ~ Tr h v\\ H ™(K) < h k+1 ~ m \v\ Hh+1{ K), m = 0,l,...,k. 

Here, in the case of a Scott-Zhang interpolation operator, K is the patch of elements 
neighboring K. 


2.3. Stability and Convergence. In this section, we prove that the finite element method 
proposed in (2.13) is stable. This is done by first proving the coercivity and continuity 
of ah defined in (2.15), followed by proving that bh defined in (2.18) satisfies the inf-sup 
condition. Combining these results proves stability of Ah. This strategy is similar to what 
can be found in na and m In particular, Verfiirth’s trick [22] is used to prove inf-sup 
stability. For a general overview of the saddle point theory used, see [3, 4[ [6]. We conclude 
the section by proving an a priori error estimate. Before we begin, we state appropriate 


norms. 
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2.3.1. Norms. In the analysis that follows, we shall use the following norms: 

i 

(2.25) \\\v\\\l = ^\\D Vi \\l ht + h\\(D v )+ v E V h , 

i =0 
1 

(2.26) MS = 5>,llo,„,. 

i =0 

(2.27) \\\(v,q)\\\l=\\\v\\\l + \\ q \\l (v, q) E V h X Q h . 


2.3.2. Interpolation estimates. Using (2.24) together with the trace inequality 
h~ l \\v\\ 2 K + h\\S7v\\ 2 K , we obtain the following interpolation estimate for v E V: 


(2.28) 


\v 


Kh.v\\\h < h k \v\ H k+i {n) . 


2 

rnA 


< 


See [in] for a proof. For the pressure p, we have 

(2.29) ||p - n h p\\h < h l+1 \v\ H i+i m . 

2.3.3. Coercivity and continuity. Establishing coercivity and continuity of is straight¬ 
forward and similar to ms- 


Lemma 4 (Coercivity of ah)- The bilinear form ah (2.15) is coercive: 
(2.30) 


v 


1 2 h<a h (v,v) V v E Vh- 


Proof. Note that the overlap term ay,.o provides the control 

i 


(2.31) 

(2.32) 

(2.33) 


I] \\ D Al Ki = ll^o||n M \ ni + IPvolln^nnx + ll^illn 


i=0 


1 

S lpvo||n M \n 1 + \\ D ( v o - t ’i)||n ho nn 1 + II^MlIn 
i 

< W Dv Wk +a hP (V,V), 
i=0 


where we have used that Qh,o \ Qi — Qo and 0 9) C hi as described in the section 
on the composite mesh above. We also note that for each element K that intersects an 
interface segment T we have the inverse bound 

(2.34) h\\(Dv) ■ n\\ 2 KnT < \\Dv\\ 2 k 

independent of the particular position of the intersection between K and T (see [TO]). 
Combining these two estimates with the standard approach to establish coercivity of 
Nitsche method (see for example [TO]) immediately gives the desired estimate. 


Lemma 5 (Continuity of ah). The bilinear form ah (2.15) is continuous: 
(2.35) a h (v,w) <\\\v\\\ h \\\w\\\ h Vv,wEV h . 


Proof. A proof in absence of a^o is found in [10]. Bounding ah,o is straightforward using 
the Cauchy-Schwarz inequality and the fact that ||Diy|| < |||iu|IU f° r an y □ 


□ p 
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2.3.4. Stability. Showing stability of the proposed finite element method involves several 
steps. First we show a preliminary stability estimate for A^ (2.13). Then the so called 
small inf-sup condition for bh (2.18) is shown using a decomposition of the pressure space 
into L 2 orthogonal components. For each of these components we show that an inf-sup 
condition holds. This is then used to show the big inf-sup condition for Ah- 

Lemma 6 (Preliminary stability estimate for Ah). It holds 


\u\ 


\i+h 2 \\mi hA „ h , 0 


< 


\u\ 


+ / i 2 ||Aw-Vp ||^ oWj0 


< 


A h ((u,p),(u,-p)). 


(2.36) 

(2.37) 

Proof. Recall the inverse estimate 

(2.38) A Ch m l \\v\\ H m^ 

(see [5], Section 4.5) where K e JCh,t and v € Vh- The first estimate in the lemma follows 
by adding and subtracting A u, using the triangle inequality and (2.38) as follows: 

h 2 \\\/p\\ 2 K <h 2 \\Vp-Au\\ 2 K + h 2 \\Au\\ 2 K 


(2.39) 

(2.40) 


<h 2 \\Vp-Au\\ 2 K +\\Du\\ 2 K , 


for each element K e /C^. The second estimate follows immediately using coercivity (2.30) 
since 


(2.41) 

(2.42) 

(2.43) 


A h ((u,p), (u, -p)) = a h (u, u) + d h ((u,p), (u, -p)) 

= a h (u, u) + h 2 \\S7p - A , u||n h]0 \ Whi0 


^ \\\u\\\l + h 2 \\Vp-Au\\ 2 


o\ u h,0 ‘ 


□ 


The pressure space can be written as the following L 2 -orthogonal decomposition: 

(2.44) Q = Q c © Qq © Qij 

where Q c is the space of piecewise constant functions on the partition (f2 ? }^ =0 of fl with 
average zero over hi and Qi is the space of L 2 functions with average zero over flj. We next 
show inf-sup conditions for Q c and Qq. Recall that the inf-sup condition for Qx is already 
established by Assumption |B| 

Lemma 7 (Inf-sup for Q c ). For each q e Q c there exists a w c 6 R with [||iw c [||/i = ||g||h 
such that 

(2-45) \\q\\ 2 h <b h (w c ,q), 

where the bound is uniform w.r.t. q. 


Proof. We first note that Q c is a one-dimensional vector space spanned by 


X = 


l^ol " 1 

-Iflil' 1 


in O 0 , 
in fR. 


(2.46) 
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Figure 5. The ball B R {x 0 ) cTflF 


Second, we note that since h2 0 an d are nonempty, there exists a nonempty open set 
U C Q such that T D U 7 ^ 0. Let now *0 € T D U be a point on the interface T and let 
B r (x 0 ) be a ball of radius R centered at x 0 as in Figure |5} The radius R is chosen such 
that B R {x 0 ) C U independently of the mesh size h. 

Now, let 7 = r nBjj(io) and note that on 7 , both the interface normal n and the jump 
[y] are constant. (In fact, [y] = Xi ~ Xo = — (|f)i | _1 + |0 2 1is constant on the entire 
interface T.) 

To construct the test function w c G Q c , we now let ip be a smooth nonnegative function 
compactly supported on B R (x 0 ) and take v(x) = cip(x)n[x], where again we note that 
both n and y are constants. The constant c is chosen such that 

(2-47) ((n ■ v), [x]) 7 = — ||y||ri- 

Integrating by parts and noting that y is constant on each subdomain 0^, i — 0,1, it 
follows that this construction of v leads to the identity 

(2-48) b h {v, y) = — «n ■ v), [y ]) 7 = ||y||£. 

Now, let w = 7 ThV G Vft,. It follows that 


(2.49) 

(2.50) 

(2.51) 

(2.52) 

(2.53) 


h{w, y) = b h {v, y) + b h (w - v, y) 

= llxlln - ((n- (w ~ v)), [y]) 7 
= llxlln — c((7ThP — p>), [x] 2 ) 7 
>Wx\\ 2 n-cCh\\x\\l 

> Ml 


The last inequality holds for all h G (0, h\ with h sufficiently small. (Note that the constants 
c and C do not depend on q.) The first inequality follows by noting that 


(2.54) 

(2.55) 

(2.56) 

(2.57) 


l((w - <p), [x] s 


< ||(7r /l ( / 9-(^)|| 7 ||[y] 2 || 7 

IIW - - <p\\m(B R ( Xo) ) 

h{\\Vip\\ B r{x 0 ) + II ll^ji(sco)) 

<hMl 


< 

< 

r^j 
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Here we have used the Cauchy-Schwarz inequality, a trace inequ ality on B R (x 0 ), an in¬ 
equality of the type ab < a 2 + b 2 , the interpolation estimate (2.24) and the definitions of x 

21 


(2.46) and (p. We also note that the estimate ||[x] J || 7 < ||x||o follows since x is (piecewise) 
constant. 

Finally, since q G Q c , we may write q = C\X for some c\ > 0. (If c\ < 0, we may redefine 
X-) Taking w c = c 2 w, where c 2 > 0 is chosen such that |||u; c |||/i = ||(?||h, we have 


(2.58) 

(2.59) 

(2.60) 

(2.61) 

since Ci = ||g|U/||x||o and c 2 = 


b h (w c ,q ) = dc 2 b h (w,x) 


^ CiC 2 

= % 
Cl 


> 


12 

I h") 


MIU and thus Ci/c 2 = |||w|||ft/||x||n ~ 1- 


□ 


Lemma 8 (Inf-sup for Q 0 ). For each g 6 Q 0 there exists a w e Wh ,o C Vfi t o with 
\\Dw\\ Uh0 = ||g - A^ 0 (g)||^ 0 such that 


2 


h 2 \\Vq\\l hA „ hfi <b h {w,q), 


we 


(2.62) ID - A n „( g )|| 

where the bound is uniform w.r.t. q. 

Proof. Recall the definitions of W^ i0 and A^ 0 from Assumption |Bj We first show that 
can change the average from Aq 0 (q) to A Wft 0 (q) using the following estimates: 

(2.63) || Q An 0 ((g)||n h0 ^ \\q ^u>h,o (?)Il^i/i.o T IIAw hi0 ( ( g) Ao 0 (<g)||o hj0 

(2.64) = Ik - -G.cMIk.o + l|Ano(A„ u (?) - ?)||n», 0 

(2.65) <lk-A™ >0 (?)|k, 0 , 

where we first added and subtracted Xu, h0 (q) and used the triangle inequality, then used 
the identity A u h0 (q) = An 0 (A Uh0 (q)), which holds since Xq 0 is an average, and finally we 
used the L 2 (Qh,o) stability |An 0 (v)| < ||u||n h0 of the average operator. 

Next we have the estimate 

(2.66) Ikllh. < Ikllho + h2 l|V9llh,V„ v 9 e «0, 

which follows by first observing that this inverse inequality holds: 

(2-67) lkll?q<^l|Vg||^+/i||g||| 12 

( 2 - 68 ) <h 2 \\^q\\ 2 Kl + \\q\\ 2 K2 , 

where K\ and K 2 are two neighboring elements sharing the face F\ 2 . Then, starting with 
w“ 0 = uih, o, we define a sequence of sets u>% 0 , n — 1 , 2 ,... consisting of the union of q 1 

and all elements K C Qh,o \ UJ h o 1 that share a face with an element in cu^q 1 . It then follows 

from (2.68) that 

(2.69) 


U <\\q\\y^+h 2 \\Vq\\i n v „-4, 

M u h,0 


n = 1,2,... 
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Using the assumption that 0 is close to fl 0 (2.12) together with shape regularity and 
quasi-uniformity of the mesh we conclude that uj% 0 = Tlh,o for some n < C for all h £ (0, h] 


where the constant is independent of h. Now (2.66) follows from a uniformly bounded 


number of iterations of (2.69). 


Combining (2.65) with (2.66), we obtain 


(2.70) 

(2.71) 

(2.72) 


Ik - 


,(9)lln M ^ll9-A 




< 


< 


Ik-A 




2 

U h,o 
2 
^>h 


, 0 + A 2 ||Vg ||^ i0V;ii0 


b h (w 0 ,q) + h 2 \\Vq\\l hAuh0 , 


where we used the fact that VA u h0 (q) = 0 and at last the inf-sup condition (2.9) to choose a 
wo e W hfi with ||Dio 0 || Wh( 


>t,o Ik A 


^h, 0 


u h ,0 such that b h {w 0 ,q ) = ||g- ^ h , 0 {q)\\l h0 - □ 


We now combine the inf-sup estimates for Q c and Q o to prove an inf-sup estimate for 

Qh- 


Lemma 9 (Small inf-sup). There are constants c > 0 and M > 0 such that for each q £ Qh 
there exists a w £ Vh with |||iu|||/i = ||g|U such that 

(2-73) Hkllh - C/i 2 ||Vg 0 ||n M \ WM - ?)> 

where the bound is uniform w.r.t. q. 

Proof. Take w c as in Lemma [7j w 0 as in Lemma [8] and W\ £ Wf h Consider the test 
function w = 5iW c + w 0 + W\ where <5i > 0 is a parameter. By writing q = q c + q 0 + q 1 £ 
Qc®Qo@Qi, we h ave 

(2.74) b h (SiW c + w 0 + w 1 ,q)= 5 1 b h (w c , q c ) + 5ib h (w c , q 0 ) + Sib h (w c , qf) 

+ b h (w 0 , q c ) +b h (w o, g 0 ) + b h (w 0 , g x ) 

=o =o 

+ b h (wi, q c ) + b h (w i, g 0 ) +b h (wi, gi) 

=o =o 

(2.75) > 5im c ||g c ||fc - 5i\b h {w c , g 0 )| - 5i|k,(ui c , gi)| 

+ m 0 |ko - An 0 (go)||o h , 0 “ Ch2 \\ Vgo||n M \ Wh>0 
+ mi||gi - An 1 (g 1 )||n 1 

(2.76) = ★. 

Note that b h (w i} q c ) = 0, i = 0,1. This follows from integration by parts since q c £ Q c , 
which is piecewise constant, and since Wi £ Wh,i > which is zero on the boundary. The 
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second term and third terms on the right-hand side can be estimated as follows 

(2.77) Si\b h (w c ,qi)\ = 5i\b h (w c ,qi - An M (g»))| 

(2-78) < 8i\\Dw c \\ nh J qi - An M (gi)||n M 

(2-79) < <5i||9c||n hi< ||9i - Af li (5 i )||n hii 

(2-80) < SfS^WqcWl + S 2 \\qi ~ A n .(ft 


H)\\n hA i 


where i — 0,1 and h 2 > 0 is a parameter. Here we have used the bound || diw|| < ||-Dv||, 
the definition of w c from Lemma 7] and the inequality ab < ea 2 + (4e) _1 6 2 , which holds for 
any e > 0. Continuing from (2.76), we use (2.80) to obtain 


(2.81) 


(2.82) 


★ > hi (m c - C5 ih 2 x ) \\q c \\ 2 h + ^(m, - C5 2 )\\qi ~ An M (g. 
- Ch 2 \\X7q 0 \\n KO \ Uhfi 


i) nn fc 


i =0 


> 


m ( hc\\l + y - A ni(9)lln hii ) - ^ 2 ||Vg 0 ||^ oWo , 


i =0 


= IH \h 

where we first choose h 2 sufficiently small and then hi sufficiently small to ensure that the 
two first terms are positive. 

Finally, we note that by construction 


(2.83) 


(2.84) 

(2.85) 


kill 1 ;$ llkclllh + ^lll^lllfc 

1=0 

1 

= hc\\l + 22 Ilk 


12 

I ^ h,i 


2—0 


and thus ||k||U ~ IklU- The desired result now follows by setting w = ||g||/ l ('iu/|||'iu|||/ l ), 
which gives 


( 2 . 86 ) 

(2.87) 


b h (w,q) = 


w \\\ h 

£ h(w,q) 


bh(w,q ) 


□ 


Proposition 10. (Big inf-sup) It holds 


( 2 . 88 ) 


(u, 


< 


sup 

(v,q)&V h xQ h 


Ah((u,p), (v,q)) 


{V, 
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Proof. Given p e Q^, take w & Vh be as in Lemma [9} First note that for we have the 
estimate 


(2.89) 

(2.90) 

(2.91) 

(2.92) 

(2.93) 


|4((u,p),(«’,o))| 

^ Vp 11 n M W,o 11 11 W, 0 

~ ^ (ll^*^l|f2/i,o\ w h,0 ll^7 , llw,oV‘ J h,,o) II II £lhfi\ u h,0 

< (ll^ll ^hfi\ u h,0 + h\\Vp\\ n 

h,o\^h,0 ) ll'^'^IIU/ li o\^Ji,0 

< 4" 1 (|pti||^ ioWo + h 2 ||Vp||^ oVM ) + 4||^||n MVM 

lft + ^ 2 ||Vp||n M v M ) + 4IHIL 


< 4 

r\^ A 


-i 


it 


where we have used the Cauchy-Schwarz inequality, the triangle inequality, the inverse 
estimate (2.38), the definition of the energy norm (2.25) and the definition of w in Lemma 

El 

Next for 4 > 0 we have 

(2.94) A h ((u,p),(u,-p) + 5 i(tw,0)) 

= A h ((u,p),(u, -p)) 


(2.95) 


+ 5 1 (a h (u,w) + b h (w,p ) + d h ((u,p ), (iy, 0)) 

> \\H\\l + h 2 \\wp\\l hoXuJho 
-4 (^IIMIIh + ^lli® 


4 (Ibllh “ ^ l|Vp||n hi0 \ WM 


-8,5. 


-i 


1°2 


\U\ 


lJ[ + ^l|Vp||^ oW , 0 -sMpU 


(2.96) 


> 

r^j 


(1 - C6^') |||u|||l + (1 - CS&') fe 2 ||Vp||^ oWo 

+ h (1 - cs 2 ) IIpIII, 


where we have used Lemmas Dill as well as ( |2.93[ ). Choosing hrst 52 sufficiently small 
and then 5i sufficiently small, we arrive at the estimate 


(2.97) 

(2.98) 

We now note that 


(«,p)IIIM||«|||£ + M 


< 


A h {{u,p),(u, -p) + 5i(tn,0)) 


(2.99) 

( 2 . 100 ) 
( 2 . 101 ) 
( 2 . 102 ) 


lll( w + Siw, -p)\\\l = |||w + 5ity|||h + Ml 

< IHIIh + 4IMl|j[ + Ibllfc 
<\M\1 + \\p\\1 

= III(«,p)IIII 
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and thus the desired estimate (2.88) follows since 

A h {{u,p),{u + S 1 w, -p))\\\h 


(2.103) 

(2.104) 


{u,p)\\\ h < 
< 


IIIMIIU 

Ah((u,p), (u + 5iw, -p)) 

|||(« + 5iw, 


□ 


2.3.5. A priori error estimate. In this section we use the approximation properties of the 
finite element spaces to show that the proposed method is optimal. 

Theorem 11. It holds 

(2.105) |||(u,p) - {u h ,Ph)\\\h < h k {\\u\\ H k+i {n) + |b||^(n)). 

Proof. By the triangle inequality we have 

(2.106) |||(w,p) - (u h ,p h )\\\ h < |||(u,p) - (7VhU,7T h p)\\\ h + \\\(n h u,TT h p) - {u h ,p h )\\\ h . 


From the approximation property o(2.28), we obtain an optimal estimate of the first term. 


To show an optimal estimate for the second term we recall the big inf-sup estimate Propo¬ 
sition m 

(2.107) \\\(n h u,TT h p) - (u h ,p h )\\\ h < A h ((n h u,n h p) - ( u h ,p h ), (v h ,q h )) 

(2.108) = ^((77^14, 7 v h p) - ( u,p ), (v h , q h )), 

where we have used a pair ( Vh,qh ) such that |||(vh) 9/i)IIU ~ 1 i 11 H ie inequality and the 


Galerkin orthogonality (2.21) to obtain the equality. The terms in Ah (2.14) may now 


be estimated individually. The optimal estimate for ah (2.15) follows immediately from 


continuity (2.35). For bh(u — nhU,qh) (2.18) we have 


(2.109) 

\b h {u - 7v h u,q h )\ 


< (^ || div(u-77/^)11^11^||^, 
\ 2—0 

(2.110) 

< (\\D(u — 77hw)||n oU n 1 ||?h||n oU f2 1 

(2.1H) 

< (|||«-77 /l tt|||^||g ft ||^ + /l|||tl-7 


1/2 


|2 \ 1/2 


uilk/illr) j 

where we have used the Cauchy-Schwarz inequality in the first two inequalities and the 


definition of the energy norm (2.25) in the last inequality. Using a similar argument we 


obtain the following estimate for bh{vh,p — 7r hp): 

(2.112) | b h {v h ,p- TT h p)\ < lllv/JHfcllp - lT h p\\ h 












16 


AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG 


(see mi)- Finally, we estimate dh to obtain 


(2.113) 


(2.114) 

(2.115) 

(2.116) 


dh({u - TT h U,p-7T h p), (v h ,q h )) 


< h 2 (J|A(u - 7r fc u)||£ MWi0 + ||V(P - n hP)Wci hy0 \uj h0 

1/2 


1/2 


< 


< 


x ^IIAv/illn^oV^o + ll^lln h>0 W,o 

\\ D ( U - 7t /i u)||^ mVm + ||p - n hP \\l h ^ ht 
2 2 \ 


1/2 


X 


U 


7r h«lllh+ IIP-’TftPlIh) 172 (IIMI/i + Uplift) 


|2\ 1/2 


= [ U 


KhulWl + Up - TT h p\\l) 1/2 |||(vh,g fc )|||ft, 


where we have used the Cauchy-Schwarz inequality, the triangle inequality, the inverse 
estimate (2.38), the definition of the energy norm (2.25) and at last the definition of the 
full triple norm (2.27). The a priori estimate now follows from the interpolation estimates 
(12224} and ((27281) □ 


3. Results and discussion 

3.1. Numerical results. To illustrate the proposed method, we here present convergence 
tests in 2D and 3D as well as a more challenging problem simulating flow around a 3D 
propeller. The numerical results are performed using FEniCS [Tl, 15], which is a collection 
of free software for automated, efficient solution of differential equations. The algorithms 
used in this work are implemented as part of the “multimesh” functionality present in the 
development version of FEniCS and will be part of the upcoming release of FEniCS 1.6 in 
2015. 


3.1.1. Convergence test. As a first test case, we consider Stokes flow in the domain = 
[0, l] d , d — 2, 3, with homogeneous Dirichlet boundary conditions for the velocity (no-slip) 
on the boundary. For d — 2, the exact solution is given by 

(3.1) u(x,y ) = 27rsin(7nr) sin ( 717 /) • (cos(7n/) sin(7rx), — cos(7rx) sin(7rp)), 

(3.2) p(x,y) = sin(27nr) sin(27rp), 


with corresponding right-hand side 


(3.3) 


f sin(27n/)(cos(27nc) — 27t 2 cos(27 tx) + 7t 2 )\ 
\ v sin(27nr)(cos(27rp) + 27 t 2 cos(27n/) — vr 2 ).y 


For d — 3, the exact solution is 


(3.4) u(x, y, z ) = sin(7rp) sin(7rz) • (1, — sin(7n/) cos(7r^), sin(7rz) cos(7rp)), 

(3.5) p(x,y,z ) = 7rcos(7rx), 
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Figure 6. Location of the overlapping domain in the background mesh. 
The domain Qj is placed in the center of D and rotated along the z-axis in 
2D (left) and along the y- and z-axes in 3D (right). 


with corresponding right-hand side 

( sin(7ry) sin(7rz) — sin(7nr)\ 
sin(27r^)(2 cos(27ry) — 1) J . 
sin(27n/)(l — 2 cos( 27 t^) J 

In both cases, the velocity field is divergence free and the right-hand side has been chosen 
to match the given exact solutions. We let the overlapping domain Di be a d —dimensional 
cube centered in the center of D with side length 0.246246 rotated 37° along the z-axis. 
For d = 3, Oi is rotated the same angle along the y-axis as well. The domains Dj are 
illustrated in Figure [6j 

The discrete spaces are P^-Pk-i Taylor-Hood finite element spaces with continuous 
piecewise vector-valued polynomials of degree k discretizing the velocity and discontinuous 
scalar polynomials of degree l = k — 1 discretizing the pressure. These spaces are inf-sup 
stable on the uncut elements of the background mesh discretizing fl 0 and on the whole of 
fli and therefore satisfy Assumption |B| 

Figures [ 7 ] and [ 8 ] show the convergence of the error in the II ( \ - and L^-norms in 2D and 3D 
respectively. Optimal order of convergence is obtained, although limited computer memory 
resources prevented a study for higher degrees than k — 3 in 3D. In the convergence plots, 
results for small mesh sizes, roughly corresponding to errors below 10 - ' have been removed 
because errors could not be reliably estimated due to numerical round-off errors in the 
numerical integration close to the cut cell boundary. 

3.1.2. Flow around a propeller. To illustrate the method on a complex geometry we create 
a propeller using the CSG tools of the FEniCS component mshr PI see Figure II (top 
left). The lengths of the blades are approximately 0.5. Then we construct a mesh of the 
domain outside the propeller, but inside the unit sphere. This is illustrated in Figure [9] 
(top right). The mesh is constructed using TetGen [21] and is body-fitted to the propeller. 
To simulate the flow around the propeller, the mesh is placed in a background mesh of 
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Figure 7. Convergence results, 2D. A rotated square is embedded in the 
unit square background mesh. Results in L 2 (left) and Hq (right) norms. 


L 1 norm convergence in 3D 



• 


/ 

/ P 

* 

/ 

M 

I* 




• • k=2, p =3.09 
■ ■ k=3, p=4.05 


10 ° 


Hq norm convergence in 3D 



Of 

*' 

* 

/ 

' P 

m 

■ 

* 

'' 



• • k=2, p = 1.98 
■ ■ k=3, p=3.06 


Figure 8 . Convergence results, 3D. A rotated cube is embedded in the unit 
cube background mesh. Results in L 2 (left) and Hq (right) norms. 


dimensions [—2, 2] 3 , where we have removed the elements with all nodes inside a sphere of 
radius 0.9, see Figure [9] (bottom). 

The simulation is setup with the inflow condition u(x, y, z ) = (0, 0, sin(7r(a;+2)/4) sin(7r(i/+ 
2)/4)) at z = —2, the outflow condition p = 0 at z = 2 and u(x,y,z) = 0 on all other 
boundaries, including the boundary of the propeller. The resulting velocity held using 
degree k — 2 is shown in Figure 10 Note the continuity of the streamlines of the velocity 


going from the finite element space defined on the background mesh to the finite element 
space defined on the overlapping mesh surrounding the propeller. 


4. Conclusions 

The finite element formulation for discretization of the Stokes problem presented has 
been demonstrated to have optimal order convergence, first by an a priori error estimates 
and then confirmed by numerical results. The finite element formulation studied in this 
work allows inf-sup stable spaces for the Stokes problem to be stitched together from 































































HIGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 


19 



FIGURE 9. Propeller geometry and meshes. Propeller geometry (top left) 
and body-fitted mesh (top right). Non body-fitted background mesh and 
propeller (bottom). 



Figure to. Flow around propeller. Colors indicate speed. 
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multiple non-matching and intersecting meshes to form a global inf-sup stable space. The 
method has several practical applications and one such prime example is the discretization 
of flow around complex objects. Future work includes the extension to time-dependent 
problems and to fluid-structure interaction. 
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